Abstract. A radiative hydrodynamic simulation for a 
typical, powerful high redshift radio galaxy is presented. 
The jet is injected at one third the speed of light into a 
10000 times denser, homogeneous medium. In the begin- 
ning of the simulation, the bow shock consists of a spher- 
ical shell that is similar to a spherical blast wave. This 
shell cools radiatively down to » 10 4 K, providing after 
6 x 10 6 yrs a neutral column of 3.8 x 10 21 cm -2 around 
the whole system. The shell starts to fragment and forms 
condensations. This absorbing screen will cover a smaller 
. and smaller fraction of the radio source, and therefore the 
' emission line region, and eventually form stars in typically 
, 10 4 globular clusters of 10 6 M Q . Approximately 10 9 M Q are 
5h ' entrained into the radio cocoon. This gas, cooling and illu- 
, minated from the radio source, could be the emission line 
' gas observed in high rcdshiftcd radio galaxies and radio 
I , loud quasars. The neutral column behind the bow shock 
fSJ ' can account for the absorption found in almost all of the 
, small sources. The globular cluster excess of ss 10 4 systems 
' found in present day brightest cluster galaxies (BCGs), 
■ which are believed to be the vestiges of these objects, is 
consistent with the presented scenario. 
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1. Introduction 

Collimated radio-luminous outflows are by now detected 
from the cosmological neighborhood up to redshifts in ex- 
cess of five. According to the unified scheme these jet- 
sources come in two classes: Radio-loud quasars, where 
the current redshift record is z = 5.8 (Fan et al. 2001), 
and radio galaxies. Among the latter ones the highest 
known redshift is z = 5.2 (van Breugel et al. 1999). In- 
frared observations indicate that radio galaxies contain 
a luminous quasar hidden by a dusty torus, at least at 
redshift above z w 0.8 ( Meiscnhcimcr ct al. 2001). The 



dominant energy source for ionization to the bow shock of 



maximum power of rad io galaxies increases with redshift 
( Dc Brcuck ct al. 2000 ), the suspicion being that the most 
powerful active ga laxies at z > 2 develop into BCGs 
(ICarilli et al. 2000| , and references therein). This idea is 
supported e.g. by the fact that the co-moving space den- 
sity of clusters at low red shift is com parable to that of 
powerful AGN at z = 2.5 flWest 1994| ). 

There is good reason to believe that the jets of higher 
redshift sources bore into a medium of higher density. If 
one just considers the cosmological expansion, the num- 
ber density of the intergalactic material (IGM) scales like 
^igm oc (1 + z) 3 , on large scales. This means e.g. that if 
the well known radio galaxy Cygnus A, which has a sur- 
rounding IGM number density of a few times 10~ 2 cm~ 3 
( parilli fc Barthcl 1996[ ), was located at redshift z = 3, 
its jet had to bore into a medium with n i=s 1 cm~ 3 . 
However, the environments of the most powerful radio 
galaxies, which are believed to highlight the highest den- 
sity peaks in the early universe, are not a priori expected 
to simply follow the Hubble flow. For example, it has to 
be taken into account that the associated galaxies ap- 
pear as bright ellipticals only at redshift z < 1.5, whereas 
above this redshift one finds clumpy regions, distributed 
over upto 100 kpc ( Pentericci et al. 2001 ), and a consid- 
erable fraction of the gas is transformed into stars around 
that redshift. Direct estimates from emission line obser- 
vations indicate IGM densities in the range of 0.1 cm~ 3 
(van Ojik ct al. 1997) upto 10 cm~ 3 , if one ascribes the 
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the jet (Bicknell et al. 2000). If the jets themselves consist 
in electron proton plasma and are subrelativistic on large 
scales - e.g. in Cygnus A a kpc scale jet velocity (vj) of 0.4 
c is favoured (Carilli & Barthel 1996, c being the velocity 
of light) - the number density in the jet (nj) is given by 



Mi 



A o Mtc ( kpc 
1.3 x 10" 4 cm~ 3 ' 1 



3wj 



(1) 



Here Mj is the jet's mass flux rate (in solar masses per 
year), Rj is the jet ra dius (about 0.6 kpc in Cygnus A, 
ICarilli fc Barthel 1996| ), and tor- the mass of hydrogen. 
Adopting an IGM number density of 1 cm -3 for high red- 
shift radio sources gives a density contrast r\ = wj/wigm ~ 
10 -4 . Due to the lower environment density (Daly 1995), 
7] is thought to be higher in low redshift sources. Because 
the expansion velocity of the radio structure scales pro- 
portional to y/fj, the average extension of radio galax- 
ies with z < 1.5 of 100 kpc versus 10 kpc at z > 2 
( Carilli et al. 2000[ ) could be nicely explained by such a 
cosmological IGM density increase, if the average active 
time stays constant. Another indication for high density 
environments at high redshift is the observed bending of 



the kpc scale radio jets (e.g. Pentericci et al. 1997). This 
underlines the clumpy nature of the IGM around high red- 
shift radio galaxies (HZRG). 

The optical emission of radio galaxies is often aligned 
with the jet axis. While below z = 1.5 this emis- 
sion is highly polarized, indicating mainly scattered light 
from a hidden AGN, above z — 2 less then 2% po- 
larization is found and stellar absorption lines are de- 
tected ( Carilli et al. 2000| ). This could indicate that - at 
high redshift - jets induce star formation. The aligned 
light could also reflect the large scale matter distribution 
flWest 1994[ ). 

Radio galaxies are surrounded by a Ly-a halo. The 



size of this halo grows with redshift (Carilli et al. 2000). 
Whilst sources at z < 1.5 have halos of 10 8 M© 
and an extension of 10 kpc, z > 2 sources have 
100 kpc halos of 10 9 M©. The width of the emission 



lines is typically (1000 ± 500) km/s ( |van Ojik et al. 1997 
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Dc Breuck et al. 2000). HZRG with diameters less than 
50 kpc show associated narrow absorption lines with a 
width of (« 40 ± 30) km/s flvan Ojik et al. 1997Q. The 



absorption is often saturated with column densities of 
(10 18 — 10 20 )cm~ 2 , and is preferentially blue-shifted by 
upto 250 km/s. These absorption systems were thought 
of to consist of sw 10 12 solar system sized clouds by 
van Ojik et al. (1997). Advanced modeling with pho- 
toionization codes lead to the proposition that the ab- 
s orber was a low den sity region, surrounding every HZRG 
( Binette et al. 2000 ). None of these models take into ac- 
count the hydrodynamic facet of the problem. This letter 
is therefore dedicated to an understanding of the above 
mentioned features based on hydrodynamic modeling of 
HZRGs. 



2. A hydrodynamic model for high redshift radio 
galaxies 

One of the basic differences between jets at low and high 
redshift is the environmental density. We consider a - typ- 
ical - jet with a density contrast of ?/ = 10~ 4 , a velocity 
of vj = 10 10 cm/s (internal Mach number 85), a jet ra- 
dius of 1 kpc, and a density of the homogeneous external 
medium of 1 cm -3 . This results in a kinetic jet luminos- 
ity of 10 45 4 erg/s, which is transfcrcd by some efficiency 
factor to radio luminosity. The velocity of advancement 
of the jet's head - and the bow shock in front of it - is 
expected to be vb ~ ^/rj «j = 1000 km/s. This velocity 
causes the external gas to get heated to a temperature of 
T = 1.38 x 10 7 K (u B /1000km/s) 2 , according to standard 
hydrodynamic shock jump conditions, neglecting the in- 
ternal energy of the pre-shock gas. The cooling time of 
this post-shock gas due to thermal bremsstrahlung with 
an emissivity of ebrcms = 2.1 x 10 _27 (n/cm~ 3 ) 2 y / T/K 
erg cm~ 3 s~ 4 will be: 



''cool 
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yrs, 
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which is less than the typical lifetime of a jet. Therefore, 
cooling has a major effect on the evolution of the post- 
shock gas. This gas will cool down to high densities and 
a temperature of about 10 4 K, where the cooling curve 
drops significantly ( [Sutherland fc Dopita 1993| ). We have 
solved the hydrodynamic equations for the above men- 
tioned parameters numerically with the code NIRVANA_C 



(Ziegler & Yorke 1997), upgraded due to optically thin 
cooling ( |Thiele 2000| ). NIRVANA_C is a second order ac- 
curate finite difference scheme. Cooling comes into play 
via the equation for the internal energy density {£) : 



d£ 

at 



— + V • (£v) = -p(V • v) - A, 



(3) 



where A is the global cooling function. NIRVANA_C also 
solves the ionization and recombination equations, by a 
time implicit method, in order to allow for nonequilibrium 




Fig. 1. Grayscale plot of the simulation described in the 
text. The levels show increasing number density with a 
factor of ten between the steps (one additional level is 
drawn at 10 01 cm -3 in order to highlight the shell). The 
simulation time is 6.0 mio yrs. 




Fig. 2. Sketch of the basic results of the simulation. 



processes. For the present simulation, three species were 
used, namely ionized and neutral hydrogen, and electrons. 
Below T = 10 6 K, an average cooling curv e for solar metal- 
icity was applied ([Stone fc Norman 1993| ) . The simulation 
was carried out in axisymmetry (2.5D), and the jet is re- 
solved with 10 points. With that resolution, NIRVANA.C 
marginally resolves shocks in the jet beam and gets ba- 
sic flow parameters right. Turbulence on small scales as 
well as details of the beam, especially of the head region, 
are unresolved (Krause & Camenzind 2001). The simula- 
tion was stopped when the thickness of the shell behind 
the bow shock shrinked below the resolution, which was 



clearly the case after a simulation time of 6.0 x 10 yrs. 

The results are quite unusual compared to simulations 
of higher rj: Most of the mass accumulates into an almost 
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spherical shell behind the bow shock. The maximum den- 
sity there is 47 cm -3 , and the gas has cooled to an average 
temperature of 1.7 x 10 4 K in regions with density above 
10 cm -3 , which still cover almost the whole shell surface. 
The mass accumulated in this shell is 6.19 x 10 10 M©, 
75 % of which are neutral. This gives an average neutral 
column for the shell of TVhi = 3.8 x 10 21 cm -2 , which has 
to be compared to the amount of matter displaced by the 
shell: 6.32 x 10 10 M . It follows, that 98 % of the mass are 
accumulated in the thin shell. Therefore 2 % of the mass 
or 1.4 x 10 9 Mq were entrained into the cocoon. 

The cocoon itself has a very peculiar structure. It is 
well known, that jet cocoons broaden when lowering rj 
(Norman et al. 1983). But in this simulation, the cocoon 
is not elongated along the jet beam, forming individual 
vortices. Instead, all the vortices join into one big vortex 
with a diameter of roughly 8 kpc. The vortex shedding is 
sometimes violent: while, at the time shown, only a small 
eddy is ejected at the beam head, at other times we ob- 
served upto half the jet beam splitting and forming an 
eddy. The big vortex turns around at a maximum velocity 
of rs 0.2 -0.25 vj. 

The velocity of the jet head differs considerably from 
higher rj simulations. At early times, it is about 10 times 
higher than the expectation above. Only after about two 
million years, the jet head moves, relatively constant, at 
maximum expected pace. 

Within the cocoon, the sound speed is, on average, 
20000 km/s. Therefore, pressure differences within the co- 
coon gas are rapidly communicated, with the result that 
the pressure is almost constant. The simulation shows that 
the bow shock is, at early times, perfectly spherical. It 
looks more like a supernova bubble than a jet cocoon. In- 
deed, if_on££ranpajrath£VBkocity of a spherical blast wave 
(e.g. Wiinsch fc Palous 2001 , with constant energy ejec- 

tion): , BW = 1412 (^g^) ^ ^ km/s to 

the velocity of the jet head, one finds, that the jet head 
outruns the blast wave only at times above 1.9 x 10 6 
yrs, for the parameters of our simulation. This is what 
we observe. We rewrite the above equation for the case, 
when the energy injection comes from a jet: vbw/vj = 
0.6 (Rj/R B w) 2/3 V 1/3 - The blast wave velocity depends 
on rj 1 / 3 , whereas the jet head velocity depends on ij 1 / 2 
(light jet limit). Because of that, the spherical state of the 
bow shock can be observed only at the earliest times, for 
higher rj. 



3. Comparison to observations 

Because of the large neutral column, this shell is a good 
candidate for the absorber in HZRGs. The width of an 
absorption line caused by this absorber would be approx- 
imately the sound speed within it ( Dyson et al. 1980[ ), 
which is in our simulation about 55 km/s, on average. 
This is also in remarkable agreement with observations 



(van Ojik et al. 1997). Typically, the material that is just 
heated by the shell emits Lyman a at a few percent of the 
kinetic jet luminosity. This is less then the observed Ly a 
luminosity, indicating that the bow shock is not the main 
emitter. At late times compared to tcoob it is justified to 
approximate the bow shock with an isothermal shock de- 
scription. Following Steffen et al. (1997), we derive for the 
thickness of the spherical thin shell: 



d = 27.7 



(m/m p )(r/10kpc)(T/10 4 K) 
Kw/100km/s) 2 



pc. 



(4) 



Here, m/m p denotes the mean molecular weight per pro- 
ton mass. The sound crossing time through this shell is 
about 5 x 10 5 yrs, which is much longer than the cooling 
time (about 10 4 yrs, Sutherland fc Dopita 1993] ). There- 
fore, it is possible, that the s hell breaks up into individual 
subshells ( Dyson et al. 1980 ). This could well correspond 
to the multiple associated absorption systems, often ob- 
served in HZRGs. These absorbers are preferentially blue 
shifted, a natural feature of our model. But almost never, 
a velocity above 250 km/s is observed. Approximating the 
sideways expansion of the shell by the laws for the spher- 
ical shell also for some time after the jet head begins to 
shape the shell, one can calculate 77, given the bubble size 
(> 12 kpc) at observation time. This yields 77 = 2 x 10~ 5 , 
which means that the full range of currently known ob- 
servations can be explained if the environmental density 
is » 5 cm~ 3 . rj in the above formula could be replaced by 
rjo(R/ Rj)~ p to account for a declining density distribu- 
tion. The rjo would be lower, in that case. It is interesting 
to note here, that very similar numbers have been found 
from quite different arguments for the environment of the 
radio galaxy 4C 41.17 QBicknell et al. 20001 ). 

The bubble interior looses a considerable amount of 
internal energy via synchrotron radiation, which is not 
included in the simulation. This lowers the requirement 
for the external density. Although the observed hydrogen 
mass in these systems barly exceeds 10 9 M Q , it is likely 
that they contain significantly more mass. Given that 
HZRG host the most powerful quasars of the universe, 
we should expect black hole masses of the order 1O 9 M0. 
The gas mass should be more than that. If the gas mass 
would be of the order 10 n ~ 12 MQ, as is suggested here, 
the gas to black hole mass ratio would be similar to the 
bulge to black hole mass ratio found at low redshift. In 
that case, most of the shell's mass would be - at any evo- 
lutionary state, the systems are observed at so far - in 
cold fragments just about to form stars. 

The shell is not only thermally, but also grav- 
itationally unstable. The time, when the instability 
on the shell surface occurs for the first time is 



( Wiinsch fc Palous 2001 , again for constant energy injec- 
tion rate and for the parameters of our model): t\ = 

2.4 x 10 7 ((. r / Vnr ^2|- I ,/ 1 n-4va„//^3r„/ r , I11 -3- 1 5) yrs. Gravi- 



r/kpc) 2 ( t) /10- 4 )(3uj/c) 3 (n/cm- 3 ) 5 i 

tational and thermal instability will support each other. 
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The Jeans mass in the shell is w 10 6 — 10 7 M Q . Hence, the 
shell will form stars in globular clusters of that mass. As- 
suming an efficiency of 10%, the shell would form, roughly, 
10 4 globular cluster systems. This is in remarkable agree- 
ment to what is found in nearby BCGs, like e.g. M87 



(Harris et al. 1998): BCGs show an excess of globular clus- 
ter systems of about 10 4 — 10 5 compared to non BCG ellip- 
ticals with comparable luminosity. Sometimes, two distin- 
guished globular cluster populations are observed. Harris 
et al. (1998) consider various formation scenarios in de- 
tail and conclude that a kind of galactic wind must have 
driven out a large fraction of the galaxies gas, before it 
was able to form stars. 

When the shell fragments, its covering factor decreases 
further and the absorption vanishes. This scenario could 
nicely explain, why no absorbers are observed for galaxies 
larger than 50 kpc. Afterwards the shell would become 
visible in the optical, due to the newly born stars, which 
will also increase the ionization of the remaining parts of 
the absorbing shell. A bubble shaped star forming region 
like that can be found around the radio galaxy 1243+036 
dvan Ojik et al. 1996Q . 

These associated absorption systems also seem to 
be observed in high redshift quasars (Baker 1997; 



Baker 2001). There, they have velocities of a few thousand 
km/s, consistent with the model, if one assumes that the 
line of sight in quasars is not far from the jet axis. In addi- 
tion, this part of the shell could be accelerated by photoab- 



sorption of photons from the quasar (Falle et al. 1981). 

In the simulation, we find that w 10 9 M© of shocked 
IGM are entrained into the cocoon. This number should 
not be taken too literally, because small scale turbulence 
could change that and high resolution studies are required 
to get this number accurate. Nevertheless, it is slightly 
more than the typical observed Lyman a emitting gas 
mass. Higher resolution would also be essential in order to 
determine the mixing properties and density of this gas. 
The big vortex, we observe in the cocoon is an ideal ac- 
celerator for the entrained gas. In the simulation, the gas 
is swept along, and accumulates mainly on the left-hand 
side and along and next to the jet beam. There it is spun 
up to a velocity of « 1000 km/s. We expect, that at higher 
resolution, small scale Kelvin-Helmholtz instabilities will 
entrain mass at the boundary of the vortex. This is also ob- 
served in other simulations (Krause & Camcnzind 2001). 



Given the high velocity of this propeller, it seems quite 
likely that emission line gas could be accelerated to the ob- 
served velocity of w 1000 km/s. All that should be studied 
in more detail with simulations of higher resolution. 
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